\documentclass[a4paper]{article} 

\usepackage{listings}
\lstset{tabsize=4, %
	frame=shadowbox, %把代码用带有阴影的框圈起来
	commentstyle=\color{red!50!green!50!blue!50},%浅灰色的注释
	rulesepcolor=\color{red!20!green!20!blue!20},%代码块边框为淡青色
	keywordstyle=\color{blue!90}\bfseries, %代码关键字的颜色为蓝色，粗体
	showstringspaces=false,%不显示代码字符串中间的空格标记
	stringstyle=\ttfamily, % 代码字符串的特殊格式
	keepspaces=true, %
	breakindent=22pt, %
	numbers=left,%左侧显示行号
	stepnumber=1,%
	numberstyle=\tiny, %行号字体用小号
	basicstyle=\footnotesize, %
	showspaces=false, %
	flexiblecolumns=true, %
	breaklines=true, %对过长的代码自动换行
	breakautoindent=true,%
	breakindent=4em, %
	%escapebegin=\begin{CJK*}{GBK}{hei},escapeend=\end{CJK*},
	aboveskip=1em, %代码块边框
	fontadjust,
	captionpos=t,
	framextopmargin=2pt,framexbottommargin=2pt,abovecaptionskip=-3pt,belowcaptionskip=3pt,
	xleftmargin=4em,xrightmargin=4em, % 设定listing左右的空白
	texcl=true,
	% 设定中文冲突，断行，列模式，数学环境输入，listing数字的样式
	extendedchars=false,columns=flexible,mathescape=true
	% numbersep=-1em
}
\usepackage{cite}
\usepackage{fancyhdr}
\usepackage{lipsum}
\usepackage{extarrows}
\usepackage{ctex}
\usepackage{amsfonts,amssymb} 
\usepackage{amsmath} 
\usepackage{cases}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{framed}
\usepackage{float}
\usepackage{ulem}
\usepackage{booktabs}
\usepackage{amsthm}
\usepackage{tikz}
\usepackage{pgfplots}
\usepackage{subcaption}
\usepackage{yhmath}
\pgfplotsset{width=10cm,compat=1.9}
\geometry{a4paper,left=3cm,right=3cm,top=3cm,bottom=3cm}
\renewcommand\thesection{\arabic{section}}  
\setcounter{section}{-1} 

\begin{document}
	\begin{center}
		{\heiti {\huge 项目作业报告}}
	\end{center}
	
	\begin{center}
		{\large 求数2101\ 叶陈昊\ 3210106359}
	\end{center}

	\section{程序设计思路与框架构建}
	本项目主要解决两个问题，一个是按点的样条插值（pp-Form）的实现，一个是B样条（B-Form）
	的实现.
	\subsection{pp-Form}
	按pp-Form的性质，我们定义如下的类：
	\begin{lstlisting}[language=C++]
		template<int order>
		class ppSpline
		{
		private:
			vector<double> knots;
			vector<polynomial<order>> polys;
		public:
			ppSpline(){};
			~ppSpline(){};
			ppSpline(const vector<double>& _knots, vector<polynomial<order>>& _polys):
			knots(_knots), polys(_polys){};
			
			vector<double> get_knots() const { return knots; }
			vector<polynomial<order>> get_polys() const { return polys; }		
			int get_num_of_knots() const{ return knots.size(); }    // 得到节点个数
			void print_ppForm_data(int sample_num, ostream &out);   // 打印密集散点以作图
																		// 传入小段区间上采样点个数
			void print_ppForm();    // 打印样条多项式
			void print_plot_data(int sample_num, string &name); // 给出文件名以存储密集散点
		};
	\end{lstlisting}
	它有两个成员变量，\texttt{knots}存储给出的节点$ \{x_1,\cdots, x_N\} $，
	\texttt{polys}存储在上述$N$个节点形成的$N-1$个区间上的多项式.
	
	注意到我们今后要使用
	到大量的多项式，我们定义了一个类\texttt{polynomial}，
	接受一个模板参数\texttt{int order}.
	\begin{lstlisting}[language=C++]
		template<int order>
		class polynomial
		{
		private:
			vector<double> coeff;
		public:
			...
		};
	\end{lstlisting}
	通过系数传入，从低次到高次构建多项式（如，传入一个\texttt{vector}
	$\{0,1,2\}$，得到的多项式为$x+x^2$）.
	
	接下来我们需要实现pp-Form的插值多项式，对于三次样条插值，本项目实现了三种样条：
	完全样条complete，自然样条natural，以及二阶导数值作为边界条件的样条second.
	这里只讨论三次的完全样条情形（其余类似推导）.
	\subsubsection{一维散点情形}\label{pp1dpt}
	定义如下函数，传入节点，节点对应函数值，边界条件（边界一阶导数值），返回三次样条类.
	\begin{lstlisting}[language=C++]
		// 给出一维散点,三次样条拟合,边界条件为complete或second
		ppSpline<3> solver( const vector<double> &knots,
							const vector<double> &vals,
							const boundary_condition boundary,
							const double &left_bound, const double &right_bound );
	\end{lstlisting}
	主要的思路是模仿讲义例题3.8，作出差商表（因此这里需要用到第二章编程中写过的库
	\texttt{dividedDiff.h}），利用式子(3.7)写出类似于式子(3.15)的线性方程组，
	其中系数矩阵为三对角矩阵.解这个方程的方法需要使用\texttt{C++}中的\texttt{Eigen3}包
	，使用这个包解决线性方程组的方法，写入了文件\texttt{linearEqnSolver.h, linearEqnSolver.cpp}中.
	
	求解后，我们得到了三次样条在各个节点上的二阶导数值$M_i, 1\le i\le n$
	.以区间$[x_i,x_{i+1}]$为例，我们已经有$f(x_i),f(x_{i+1}),
	M_i,M_{i+1}$，这四个信息唯一确定了一个三次多项式.我们将其写成一个函数，
	存入到\texttt{polynomial.h}中.
	\begin{lstlisting}[language=C++]
		// 给定两个节点, 节点对应函数值和二阶导数值, 求出3次多项式
		polynomial<3> cubic_0_2_val(double a, double b,
									double fa, double fb,
									double Ma, double Mb);
	\end{lstlisting}

	综上，我们求解了一维散点情形.

	\subsubsection{拟合函数情形}
	定义如下函数：
	\begin{lstlisting}[language=C++]
		// 给出插值点,函数,边界条件,三次样条拟合函数y=f(x)
		ppSpline<3> solver( const vector<double> &knots, const Function &fun,
							const boundary_condition boundary );
	\end{lstlisting}
	这里函数采用类的方式定义，见\texttt{Function.h, Function.cpp}，主要定义了
	在$x$处的函数值，一阶导数和二阶导数值.
	
	由于函数自然给出了complete或second边界条件的值，所以不需要另传参数.
	我们只需要计算好相应边界条件所需要的值，然后传入\ref{pp1dpt}中定义的函数中即可.
	以三阶完全样条为例，给出一种实现.
	\begin{lstlisting}[language=C++]
		...
		int N = knots.size();
		vector<double> vals(N);
		for(int i = 0; i < N; i++){
			vals[i] = fun.f(knots[i]);
		}
		if(boundary == complete){
			return solver(knots, vals, boundary, 
							fun.df(knots[0]), fun.df(knots[N-1]));
		}
		...
	\end{lstlisting}
	
	\subsubsection{高维散点情形}\label{usenatural}
	定义如下函数，传入散点，边界条件（边界散点关于弧长参数的梯度），返回三次样条类的向量.
	\begin{lstlisting}[language=C++]
		// 给出高维散点,三次样条拟合,边界条件为complete或second
		vector<ppSpline<3>> solver( const vector<vector<double>> &points, 
									const boundary_condition boundary,
									const vector<double> &left_bound, 
									const vector<double> &right_bound );
	\end{lstlisting}

	采用讲义算法3.72，拆解成一组“一维散点”情形问题，用\ref{pp1dpt}中给出的函数
	解决每一个“一维散点”问题，最后得到一组三次样条，存入三次样条类的向量中.
	
	注意到要传入的是“边界散点关于弧长参数的梯度”，这在已知函数时仍然可能不好确定，
	所以一般情况下我们采用的是自然样条. 关于其他边界条件，也给出了函数接口，
	用户如果已知边界条件信息可以传入，否则用自然样条更佳.
	
	\subsubsection{线性插值情形}
	对于一维情形，已知一段区间$[x_i,x_{i+1}]$上的函数值$f(x_i),f(x_{i+1})$，
	就唯一确定了一个一次多项式.求解一次多项式的函数如下.
	\begin{lstlisting}[language=C++]
		// 给定两个节点和节点对应函数值, 求出1次多项式(线性函数)
		polynomial<1> linear_0_val(double a, double b, double fa, double fb);
	\end{lstlisting}
	对于高维情形，较三次样条相比更加简单，只需要传边界函数值向量即可.
	下面给出求解一次样条的接口.
	\begin{lstlisting}[language=C++]
		// 给出插值点,函数,一次样条(分段线性插值)拟合y=f(x)
		ppSpline<1> solver(const vector<double> &knots, const Function &fun);
		// 给出一维散点,一次样条(分段线性插值)拟合
		ppSpline<1> solver(const vector<double> &knots, const vector<double> &vals);
		// 给出高维散点,一次样条(分段线性插值)拟合
		vector<ppSpline<1>> solver(const vector<vector<double>> &points);
	\end{lstlisting}

	\subsection{B-Form}
	按B-Form的性质，我们需要按如下步骤得到B样条.
	\begin{itemize}
		\item 样条基函数：给定节点$t_1,\cdots,t_N$，得到一个紧支集为$[t_1,t_N]$的样条
		基函数（它的次数$n$可以由节点数$N$唯一确定，事实上，$n = N-2$）
		\item 基：给定节点$t_1,\cdots,t_N$和样条多项式的次数$n$，生成一组样条基函数，
		它们至少在某个给定节点上是非零值
		（即，生成关于这些节点的样条空间$\mathbb{S}_n^{n-1}(t_1,\cdots,t_N)$的一组基
		$B_{2-n}^n(x),B_{3-n}^n(x),\cdots,B_N^n(x)$，
		见讲义定理3.49）.此时可以得到基的个数为$n+N-1$.
		\item B样条：给定节点$t_1,\cdots,t_N$, 样条多项式的次数$n$, 系数$a_{2-n},\cdots,a_{N}$，结合第二点得到的基，可以唯一确定一个B样条多项式：
		$S(x) = a_{2-n}B_{2-n}^n(x)+\cdots+a_NB_N^n(x).$
		
	\end{itemize}
	
	\subsubsection{样条基函数的实现}\label{basicBfun}
	定义样条基函数类\texttt{basicBfun}如下，以多项式次数\texttt{order}作为模板.
	\begin{lstlisting}[language=C++]
		template <int order>
		class basicBfun
		{
		private:
			vector<double> knots;              
			vector<polynomial<order>> polys;    
		public:
			basicBfun(){};
			// 传入节点, 构造关于这组节点的样条插值基函数
			basicBfun(vector<double>& _knots);
			~basicBfun(){};
			...
		};
	\end{lstlisting}
	注意到样条基函数是以分段区间上的多项式实现的，因此需要用\texttt{vector}容器储存.
	具体实现见\texttt{BSpline.h}.
	
	在难度上，只要求实现小于等于3次的并不难，这里本人没想到别的方法，直接暴力穷举了.
	
	\subsubsection{基的实现}\label{basicBfuns}
	定义样条空间的基\texttt{basicBfun}如下，以多项式次数\texttt{order}作为模板.
	\begin{lstlisting}[language=C++]
		template<int order>
		class basicBfuns
		{
		private:
			vector<double> knots;               // e.g. 1 2 3 4 5 6 7 8, order = 3
			vector<double> extend_knots;        // -2 -1 0 1 2 3 4 5 6 7 8 9 10 11
			vector<basicBfun<order>> polys;    // 区间:[-2,2],[-1,3],...,[7,11]
		public:
			basicBfuns(){};
			// 传入节点, 得到延伸节点和对应样条空间的一组基
			basicBfuns(const vector<double>& _knots);
			~basicBfuns(){};
			...
		};
	\end{lstlisting}
	注意到这里有两种节点，一个是给定的节点，另一个是需要延伸出去的节点.
	我们延伸节点的方式是，取给定节点$t_1,\cdots,t_N$生成区间的平均区间长度
	$l = \frac{t_N-t_1}{N-1}$，然后以$l$长度等距地延伸出去，
	在左、右两端各延伸\texttt{order}
	个节点，得到新的节点组\texttt{extend\_knots}.样条空间的基中元素的紧支集就必然
	包含在延伸节点所形成的区间中.
	
	这里的难度不高，因为得到新节点组后，只需调用第一点\ref{basicBfun}的相关
	函数即可构造.
	
	\subsubsection{B样条$S(x)$的实现}
	要想算出$S(x)$，必然要得到基和系数，这是线性代数告诉我们的.基在\ref{basicBfuns}
	中已经可以求得.接下来的主要任务是求系数.
	
	由于本项目对B样条不要求太高的任意性，故程序只提供了两种可求解的函数：
	\begin{itemize}
		\item 三次B样条，插值点就是节点（主要考虑到，如果插值点不是节点的时候，
		首先系数矩阵不是三对角矩阵，其次对边界条件的处理会更复杂，所以只实现了这个方法），
		配有complete边界条件.
		\item 二次B样条，插值点安插在节点之中，二者交错分布（主要考虑到，如果插值点
		不是这种情形，系数矩阵也会不是三对角矩阵；当然，一般选插值点肯定均匀分布为佳），
		配有边界条件：在边界节点的值相等（$S(t_1)=f(t_1),S(t_N)=f(t_N)$）.
	\end{itemize}

	我们只讨论散点情形.对于传入函数情形，处理办法同pp-Form，这里就不赘述了.
	
	函数接口如下：
	\begin{lstlisting}[language=C++]
		// 散点拟合求解
		template <int order>
		vector<double> solver(const vector<double>& knots, 
								const vector<double>& sites, const vector<double>& vals,
								const boundary_condition boundary, 
								double left_bound, double right_bound);
	\end{lstlisting}
	其中\texttt{boundary\_condition}就是上面给出的两种可求解类型的边界条件.
	
	对于这两个B样条，结合讲义定理3.57和3.58的思路，都会得到一个线性方程组，
	系数矩阵为三对角矩阵，待求解的项为$S(x)$中的系数$a_{3-n},a_{4-n},\cdots,a_{N-1}$
	，右端项（除去头尾两项）为散点的纵坐标值（它的横坐标就是插值点），右端头尾项需要
	再加一些边界条件约束.最后
	用\texttt{linearEqnSolver.h}中函数得到解，再在解的头尾加入边界系数$a_{2-n},a_N$
	，即可得到基对应的一组系数$a_{2-n},\cdots,a_{N}$，从而得到$S(x)$.
	
	\subsection{代码测试}
	
	对一维散点的三次pp-Form，ProblemA中已经全部测试了一遍，这里不进行测试.
	
	对高维散点的三次pp-Form（自然样条），ProblemE中已测试，这里不进行测试.
	
	对B-Form，我们验证定理3.57和3.58，而ProblemC中已测试，这里不进行测试.
	
	我们主要进行高维散点的线性插值和三次pp-Form（完全样条）的测试.
	为方便测试（能够确定边界散点关于弧长参数的梯度）
	，我们对$x^2+y^2=1$进行拟合. 测试代码见\texttt{testCircle.cpp}.
	
	选取圆周上从$(1,0)$开始的6个六等分点作为散点进行拟合，拟合结果如下.
	\begin{figure}[H]
		\centering
		\begin{tikzpicture}
			\begin{axis}[
				%		title=用曲线拟合平面散点,
				legend pos=outer north east,
				axis equal, 
				xlabel=$x$,ylabel=$y$,
				]			
				\addplot+ [only marks] coordinates {  
					(1,0)
					(0.5,{sqrt(3)/2})
					(-0.5,{sqrt(3)/2})
					(-1,0)
					(-0.5,{-sqrt(3)/2})
					(0.5,{-sqrt(3)/2})
				};
				\addlegendentry{散点}

				\addplot[red, mark=none] table {output/testCircle_data/cubicCplPPS.txt};  
				\addlegendentry{三次complete样条拟合}
				
				\addplot[blue, mark=none] table {output/testCircle_data/linearPPS.txt};
				\addlegendentry{线性插值拟合}
				
				\addplot[green, mark=none, domain=0:2*pi, samples=100] ({cos(deg(x))}, {sin(deg(x))});  
				\addlegendentry{$x^2+y^2=1$}

			\end{axis}
		\end{tikzpicture}
		\caption{两种样条拟合}
	\end{figure}
	
	可以看到，线性插值拟合正确，三次complete样条拟合的效果也非常不错，说明程序
	基本能够达到预期.
	
	\subsection{补充说明与待改进之处}
	\begin{itemize}
		\item 本项目对pp-Form和B-Form均进行了一定的包装，使得用户能够通过
		输入散点、边界条件，能够得到一种对应样条的返回值.
		\item 本项目没有将pp-Form和B-Form
		这两类合起来再进行包装，让用户能够通过输入自己想要的样条形式
		得到相应的结果.主要是因为二者虽然都是一些多项式段，但二者内部的成员变量是
		不一样的：pp-Form是节点和区间上的一个多项式，而B-Form存储着样条基和对应系数
		，并且有节点和插值点不一样的情形，
		故二者不好合并起来.不过要想合并，理论上可以将后者做一个多项式加法
		（在\texttt{polynomial.h}中实现），以得到每个区间段上的多项式，使之结构
		与前者相仿.在今后可能会尝试一下.
		\item 对于ProblemA-E的输出，除了ProblemA的子区间中点的最大范数和
		ProblemD的结果直接输出在终端上，其他
		数据（基本上是样本点数据，用以作图）都是以文件写入的.
		如果成功写入，终端会有输出提示.在\texttt{ppSpline.h}和
		\texttt{BSpline.h}中配备了一些打印多项式的函数，如有需要，可以使用.
		\item 用\texttt{make}编译\texttt{cpp}文件、\texttt{make report}编译
		\texttt{tex}文件的时间较长，还请耐心等待.
	\end{itemize}
	
	\newpage
	下面进入解决问题的部分.
	\section{ProblemA}
	\subsection{三次样条拟合函数}
	
	见\texttt{ProblemA.cpp}.本题分别使用了complete，second，natural
	边界条件对$f(x)=\frac{1}{1+25x^2}$进行三次样条拟合.在每种边界条件下，
	对$N=6,11,21,41,81$进行节点选取并在上面作出样条多项式.对每个样条多项式
	采样，描点并连线，得到五个样条曲线.结合原曲线$y=f(x)$，我们作出如下的图.
	
	\begin{figure}  [H]
		\centering  
		\begin{subfigure}{0.45\textwidth}  
			\includegraphics[width=\linewidth]
			{output/ProblemA_data/complete.pdf}   
			\caption{complete边界条件}  
		\end{subfigure}  
		\hfill  
		\begin{subfigure}{0.45\textwidth}  
			\includegraphics[width=\linewidth]
			{output/ProblemA_data/second.pdf}  
			\caption{second边界条件}  
		\end{subfigure}  
		\hfill  
		\begin{subfigure}{0.45\textwidth}  
			\includegraphics[width=\linewidth]
			{output/ProblemA_data/natural.pdf}  
			\caption{natural边界条件}  
		\end{subfigure}  
 
		\caption{不同边界条件拟合的三次样条}  
	\end{figure}     
	
	可以看到，这三张图基本上没太大差别，并且样条插值不会产生Runge现象.
	
	\subsection{误差估计与分析}
	对每个固定的$N$，可以得到样条在每个子区间中点处的值，再借由精确函数在这个地方的值，
	得到二者之间的差值.通过对比这些差值的绝对值，我们以其中的最大值作为最大范数评价的标准，
	得到的结果为：
	 
	\begin{table}[h]  
		\centering  
		\caption{子区间中点的最大范数}  
		\begin{tabular}{|c|c|c|c|}  
			\hline  
			$N$ & \textbf{complete} & \textbf{second} & \textbf{natural} \\  
			\hline  
			6 & 0.421705 & 0.423703 & 0.423482 \\  
			\hline  
			11 & 0.0205289 & 0.0205291 & 0.0205306 \\  
			\hline  
			21 & 0.00316894 & 0.00316894 & 0.00316894 \\  
			\hline  
			41 & 0.000275356 & 0.000275356 & 0.000275356 \\  
			\hline  
			81 & 1.609e-05 & 1.609e-05 & 1.609e-05 \\  
			\hline  
		\end{tabular}  
	\end{table}
	
	观察发现，在节点数很多时，无论何种边界的样条插值，误差都不会有太大的区别.
	
	另外，从表中可以看出，每当子区间的个数$(N-1)$扩大一倍时，最大误差$E$会下降一个数量级.
	因此，我们可以猜测，$\ln(N-1)$与$\ln E$成线性关系.接下来以完全样条为例，验证这个猜想.
	\begin{figure}[H]
		\begin{center}
			\begin{tikzpicture}    
				\begin{axis}[    
			%		title = {$\ln E$ 与 $\ln(N-1)$ 的相关性},  
					legend pos=outer north east,  
					xlabel = {$\ln(N-1)$},  
					ylabel = {$\ln E$},  
					]    
					\addplot+ [only marks] coordinates {  
						({ln(5)}, {ln(0.421705)})  
						({ln(10)}, {ln(0.0205289)})  
						({ln(20)}, {ln(0.00316894)})  
						({ln(40)}, {ln(0.000275356)})  
						({ln(80)}, {ln(1.609e-05)})  
					};    
					\addplot [
					domain=1.5:4.5, 
					samples=100, 
					color=black,
					]
					{-4*x+6}
					node{$k=-4$};
				\end{axis}    
			\end{tikzpicture} 
		\end{center}
		\caption{$\ln E$ 与 $\ln(N-1)$ 的相关性} 
	\end{figure}
	
	
	可以看到，$\ln(N-1)$与$\ln E$呈现着很强的线性相关性，且这条直线的斜率近似为$-4$.
	故我们猜测$\ln E \approx -4\ln(N-1)$. 结合$h(N-1) = 2$（$h$为均匀步长）知，
	$$ E = C(N-1)^{-4} = C'h^4. $$
	其中$C,C'$是与$N$无关的一个常数.
	此式与讲义定理3.12中误差在$j=0$时的估计一致.从而我们验证了三次样条的误差收敛阶为4.
	
	\section{ProblemB}
	在\texttt{BSpline.h}中已经给出一般节点下的三阶complete和二阶边界节点函数值相等的
	两种B样条.此题要验证定理3.57和3.58，对于整数节点来说，当然也成立.
	具体的两种函数接口见\texttt{ProblemB.h}.
	\begin{itemize}
		\item 定理3.57：传入要进行拟合的区间的左、右端点值（整数$a<b$），以及函数$f$，
		由于节点和插值点相同，只需这3个信息就可以求得cardinal \ B样条.
		\item 定理3.58：传入$a,b,f$（同上定义）以及插值点组.
	\end{itemize}
	
	\section{ProblemC}
	采用B题的头文件.具体见\texttt{ProblemC+D.cpp}.测试结果如下：
	\begin{figure}[H]
		\begin{center}
			\begin{tikzpicture}    
				\begin{axis}[      
					legend pos=outer north east,  
					xlabel = {$x$},  
					ylabel = {$y$},  
					]    
					\addplot[red, mark=none] table {output/ProblemC_data/cubic.txt};
					\addlegendentry{3次拟合曲线}
					\addplot[black, mark=none] table {output/ProblemC_data/quadratic.txt};
					\addlegendentry{2次拟合曲线}
					\addplot[blue, mark=none] table {output/ProblemC_data/f.txt};
					\addlegendentry{$f(x)=\frac{1}{1+x^2}$}
				\end{axis}    
			\end{tikzpicture} 
		\end{center}
		\caption{用cardinal\ B样条拟合函数$f(x)$} 
	\end{figure}\label{ProblemC}
	
	从图上可以直观看出，二次cardinal \ B样条没有三次cardinal\ B样条拟合的好.
	
	\section{ProblemD}
	
	由题意，在\texttt{ProblemC+D.cpp}中打印了两种cardinal\ B样条在相同点$x$处的
	插值误差$E_S(x)$，如表\ref{ESx} 所示.
	
	\begin{table}[h]
		\centering
		\caption{两种cardinal\ B样条的插值误差$E_S(x)$}
		\label{ESx}
		\begin{tabular}{|c|c|c|}
			\hline
			$x$	& \textbf{cubic} & \textbf{quadratic} \\
			\hline
			-3.5 & 0.000669568 & 0 \\
			\hline
			-3 & 0 & 0.00141838 \\
			\hline
			-0.5 & 0.0205289 & 1.11022e-16 \\
			\hline
			0 & 0 & 0.120238 \\
			\hline
			0.5 & 0.0205289 & 0 \\
			\hline
			3 & 1.38778e-17 & 0.00141838 \\
			\hline
			3.5 & 0.000669568 & 0 \\
			\hline
		\end{tabular}
	\end{table}

	从这张表中看到，有些值在计算机的误差精度下为0，这是因为二者插值点种类导致的.
	cubic插值点为整数点，所以在$x=-3,0,3$时误差为0；quadratic插值点为半整数点，
	所以在$x=-3.5,-0.5,0.5,3.5$时误差为0.直观上，位于quadratic的$x=0$处
	的插值误差比表内其他数据都大（1个数量级及以上），也从一个方面说明拟合效果不如cubic.
	
	\section{ProblemE}
	题目要求绘制如下的心形曲线：
	\begin{align}
		x^2+\left(\frac{3}{2}y-\sqrt{|x|}\right)^2=3\label{3.88}
	\end{align}

	首先，令$x=0$，解得$y=\pm\frac{2}{\sqrt{3}}\approx \pm1.1547$.
	记$A = (0,-\frac{2}{\sqrt{3}}), B=(0,\frac{2}{\sqrt{3}})$.
	不难注意到这个心形曲线沿$AB$（即$y$轴）对称，因此我们只需要考虑心形曲线
	在$y$轴右侧的部分（记作$\wideparen{AB}$）.
	我们可以先作出这个心形曲线，以便观察如何更好地选取特征点进行样条插值.
	\begin{figure}[H]
		\begin{center}
			\begin{tikzpicture}
				\begin{axis}[
				%	title=用曲线拟合平面散点,
				%	legend pos=outer north east,
					xlabel=$x$,ylabel=$y$,]
					\addplot+ [only marks] coordinates {  
						(0,{-2/sqrt(3)}) 
						% 对称点  
						(0,{2/sqrt(3)}) 
					}; 
					\node[right] at (axis cs:0,{-2/sqrt(3)}) {$A$};  
					\node[right] at (axis cs:0,{2/sqrt(3)}) {$B$};  
					
					\addplot[red, mark=none] table {output/ProblemE_data/heartcurve.txt};
				%	\addlegendentry{心形曲线}
				\end{axis}
			\end{tikzpicture}
		\end{center}
		\caption{心形曲线}
	\end{figure}

	在图中可以看到，$A$和$B$两点是曲线上的全部尖点（不光滑点），由式(\ref{3.88})不难得到弧
	$\wideparen{AB}$上的点满足的方程：
	\begin{align}
		y = \frac{2}{3}\left(\sqrt{x}\pm\sqrt{3-x^2}\right),
		 0\le x\le \sqrt{3}.\label{rteqn}
	\end{align}

	现在我们需要在弧$\wideparen{AB}$上选取离散点并进行样条拟合.
	根据\ref{usenatural}中的描述，我们不知道$A,B$点处关于弧长参数的梯度，
	所以我们采用natural边界条件来进行样条拟合.
	
%	我们以$B$点为例.在$B$的右邻域，曲线满足解析式
%	$y := f(x) = \frac{2}{3}(\sqrt{x}+\sqrt{3-x^2})$，从而$f\in C([0,\sqrt{3}])
%	\cap C^2((0,\sqrt{3}))$.进一步计算得
%	\begin{align}
%		f'(x)=\frac{1}{3}\left(\frac{1}{\sqrt{x}}-\frac{2x}{\sqrt{3-x^2}}\right),\\
%		f''(x)=-\frac{1}{6}x^{-\frac{3}{2}}-2(3-x^2)^{-\frac{3}{2}}.
%	\end{align}
%	从而在$B$点处，$y$方向的导数和二阶导数均为$\infty$.上述对$A$点也同理.
%	故complete边界条件和second边界条件在
%	理论上不太适用.因此，我们考虑采用natural边界条件来拟合弧
%	$\wideparen{AB}$.
	
	鉴于曲线的对称性和已知的起终点$A,B$（作为两个特征点），
	我们把题目中的$n=10,40,160$约化为$n'=6,21,81(\mathrm{i.e.},n'=\frac{n}{2}+1)$，
	即只需要对弧$\wideparen{AB}$上的$n'$个散点进行样条拟合.
	
	对$n'=6$的情形，我们需要在$y$轴右侧再找四个特征点.注意到$A,B$点具有奇性，我们
	应该会选取多一些靠近它们的点.最终选定的四个特征点如下：
	$$ (0.5,-0.6341),
	(1.3028,0),
	(1.7218,1) ,
	(1,1.6095). $$
	借助这六个点（加上对称的四个点）做出如下图的拟合.

	\begin{figure}[H]
		\centering
		\begin{tikzpicture}
			\begin{axis}[
			%		title=用曲线拟合平面散点,
					legend pos=outer north east,
					xlabel=$x$,ylabel=$y$,
					%			axis lines=center,
					%			tick align=inside,
					%			xmin=-1.5,xmax=4.9,ymin=-3.3,ymax=3.9,
					%			xtick distance=1,ytick distance=1,
					%			axis equal image
					]
					\addplot+ [only marks] coordinates {  
							(0,-1.1547)  
							(0.5,-0.634137)
							(1.30277561,0)
							(1.7218382,1) 
							(1,1.6094757) 
							(0,1.1547)  
							% 对称点  
							(-1,1.6094757)
							(-1.7218382,1)
							(-1.30277561,0)
							(-0.5,-0.634137)
							(0,-1.1547)
						}; 
					
					\node[right] at (axis cs:0,{-2/sqrt(3)}) {$A$};  
					\node[right] at (axis cs:0,{2/sqrt(3)}) {$B$}; 
					\node[right] at (axis cs:0.5,-0.634137) {$C$}; 
					\node[above] at (axis cs:0.5,1.576946) {$D$}; 
					\node[right] at (axis cs:{sqrt(3)},{2/3*sqrt(sqrt(3))}) {$E$};
					\addlegendentry{节点（$n=10$）}
					%			\addplot+ coordinates {(2,1)};
					%			\addplot+ coordinates {(3,-2)};
					\addplot[red, mark=none] table {output/ProblemE_data/10knots.txt};  
					\addlegendentry{拟合曲线}
					
					\addplot[blue, mark=none] table {output/ProblemE_data/heartcurve.txt};
					\addlegendentry{心形曲线}
					\addplot+ [only marks] coordinates {  
						(0.5,1.576946)
						({sqrt(3)},{2/3*sqrt(sqrt(3))})
					}; 
				\end{axis}
		\end{tikzpicture}
		\caption{$n = 10$的样条拟合}\label{n'=6}
	\end{figure}

	当然，由于只有4个插值点可以使用，我们不能兼顾到所有情况，
	如在$A,B$点附近实际上并没有提供新的插值点，
	因为要事先保证$\wideparen{AB}$中间段的光滑性.
	在之后节点增多时，会考虑补充这一情形.
	
	对$n'=21$的情形，我们可以使用19个插值点，我们的初步想法是，
	分别在$A,B$点附近使用6个插值点，
	而在剩下的曲线中较均匀地撒上7个插值点.
	
	我们在图\ref{n'=6}中选取了如下特征点：
	$$C(0.5,-0.6341),D(0.5, 1.5769)
	,E(\sqrt{3},\frac{2}{3}\sqrt[4]{3}).$$
	其中$E$为心形曲线的最右端点.
	这样，$\wideparen{AB}$被划分成了4段小弧.
	在每段小弧上（不含端点）较均匀地撒4个点，从而我们得到了这21个插值点.
%		$\wideparen{AC}$和$\wideparen{BD}$上（不含端点）
%		较均匀地撒4个点，在$\wideparen{CE}$和
%		$\wideparen{CD}$上较均匀地撒4个点.
	最终经过这些插值点得到的拟合曲线如下.
	
	\begin{figure}[H]
		\centering
		\begin{tikzpicture}
			\begin{axis}[
				%		title=用曲线拟合平面散点,
				legend pos=outer north east,
				xlabel=$x$,ylabel=$y$,
				%			axis lines=center,
				%			tick align=inside,
				%			xmin=-1.5,xmax=4.9,ymin=-3.3,ymax=3.9,
				%			xtick distance=1,ytick distance=1,
				%			axis equal image
				]			
				\addplot+ [only marks] coordinates {  
					(0,{-2/sqrt(3)}) 
					(0,{2/sqrt(3)})
					(0.5,-0.634137) 
					(0.5,1.576946) 
					({sqrt(3)},{2/3*sqrt(sqrt(3))})
				};
				\addlegendentry{特征点}
				\node[right] at (axis cs:0,{-2/sqrt(3)}) {$A$};  
				\node[right] at (axis cs:0,{2/sqrt(3)}) {$B$}; 
				\node[right] at (axis cs:0.5,-0.634137) {$C$}; 
				\node[above] at (axis cs:0.5,1.576946) {$D$}; 
				\node[right] at (axis cs:{sqrt(3)},{2/3*sqrt(sqrt(3))}) {$E$};
				\node[right] at (axis cs:0.25, -0.809246) {$C'$};
				\node[right] at (axis cs:1.5, 0.239146) {$F$};
				\node[above right] at (axis cs:1, 1.60948) {$G$};
%					\addlegendentry{散点（$n=40$）}
				%			\addplot+ coordinates {(2,1)};
				%			\addplot+ coordinates {(3,-2)};
				\addplot[red, mark=none] table {output/ProblemE_data/40knots.txt};  
				\addlegendentry{拟合曲线}
				
				\addplot[blue, mark=none] table {output/ProblemE_data/heartcurve.txt};
				\addlegendentry{心形曲线}
				\addplot+ [only marks] coordinates {  
					(0.25,-0.809276)
					(1.5, 0.239146)
					(1,1.60948)
				};
			\end{axis}
		\end{tikzpicture}
		\caption{$n = 40$的样条拟合}\label{n'=21}
	\end{figure}

	从中可以看到，相比于$n=10$而言，$n=40$拟合的更加完美了.不过美中不足的是，
	反而光滑段$\wideparen{DE}$的拟合效果不是很好，对$A,B$点附近的拟合
	也不是很足够.我们在$n=160$时再进行补充考虑.
	
	对$n'=81$的情形，我们考虑在曲线上点的切线与$y$轴平行的点（即点$A,B,E$），
	在它们附近撒相对多的点，而在其余地方仍然稍均匀地撒插值点.
	
	我们在图\ref{n'=21}中重新选取7个特征点:
	$$A,B,C'(0.25,-0.809276),D,E,F(1.5,0.239146),G(1,1.60948).$$
	在每段弧上较平均地撒点，撒点情况如下：（均不包含弧的端点）
	
	\begin{table}[h]  
		\centering  
		\caption{曲线段上撒点情况}  
		\begin{tabular}{|c|c|}  
			\hline  
			弧 & 撒点个数 \\  
			\hline  
			$\wideparen{AC'}$ & 13\\  
			\hline  
			$\wideparen{C'F}$ & 11\\  
			\hline
			$\wideparen{FE}$ & 13\\  
			\hline  
			$\wideparen{EG}$ & 13\\ 
			\hline  
			$\wideparen{GD}$ & 11\\   
			\hline 
			$\wideparen{DB}$ & 13\\    
			\hline   
		\end{tabular}  
	\end{table}
	
	对这81个插值点做样条拟合，得到拟合曲线如下：
	
	\begin{figure}[H]
		\centering
		\begin{tikzpicture}
			\begin{axis}[
				%		title=用曲线拟合平面散点,
				legend pos=outer north east,
				xlabel=$x$,ylabel=$y$,
				%			axis lines=center,
				%			tick align=inside,
				%			xmin=-1.5,xmax=4.9,ymin=-3.3,ymax=3.9,
				%			xtick distance=1,ytick distance=1,
				%			axis equal image
				]			
				\addplot+ [only marks] coordinates {  
					(0,{-2/sqrt(3)}) 
					(0,{2/sqrt(3)})
					(0.5,1.576946) 
					({sqrt(3)},{2/3*sqrt(sqrt(3))})
					(0.25,-0.809276)
					(1.5, 0.239146)
					(1,1.60948)
				};
				\addlegendentry{特征点}
				\node[right] at (axis cs:0,{-2/sqrt(3)}) {$A$};  
				\node[right] at (axis cs:0,{2/sqrt(3)}) {$B$}; 
				\node[above] at (axis cs:0.5,1.576946) {$D$}; 
				\node[right] at (axis cs:{sqrt(3)},{2/3*sqrt(sqrt(3))}) {$E$};
				\node[right] at (axis cs:0.25, -0.809246) {$C'$};
				\node[right] at (axis cs:1.5, 0.239146) {$F$};
				\node[above right] at (axis cs:1, 1.60948) {$G$};
				\addplot[red, mark=none] table {output/ProblemE_data/160knots.txt};  
				\addlegendentry{拟合曲线}
				
				\addplot[blue, mark=none] table {output/ProblemE_data/heartcurve.txt};
				\addlegendentry{心形曲线}
			\end{axis}
		\end{tikzpicture}
		\caption{$n = 160$的样条拟合}\label{n'=81}
	\end{figure}
	
	可以看到，在对$n'=21$的情形再加以改进后，$n'=81$的情形已经达到一个几乎完美拟合
	的状态.当然如果仔细地看，会发现$A,B,E$点附近还是可以做到更加好的一个拟合效果.这里没做到
	的主要原因是，之前提到的所谓“较平均的撒点方式”，只是相对于横坐标$x$平均，所以当
	函数很陡峭时（导数达到$\infty$时），撒点就不是很平均了.理论上最好的方式当然是根据
	弧长取平均，不过本人目前没有想到很好的方法，使得我们能够简单地得到以弧长平均作为指标
	的“平均撒点”办法.
	
	最后，我们结合$n=10,40,160$，以及原心形曲线，将它们绘制在一张图上.
		
	\begin{figure}[H]
		\centering
		\begin{tikzpicture}
			\begin{axis}[
				%		title=用曲线拟合平面散点,
				legend pos=outer north east,
				xlabel=$x$,ylabel=$y$,
				]
				\addplot[red, mark=none] table {output/ProblemE_data/10knots.txt};  
				\addlegendentry{$n=10$}
				
				\addplot[green, mark=none] table {output/ProblemE_data/40knots.txt};  
				\addlegendentry{$n=40$}
				
				\addplot[brown, mark=none] table {output/ProblemE_data/160knots.txt};  
				\addlegendentry{$n=160$}
				
				\addplot[blue, mark=none] table {output/ProblemE_data/heartcurve.txt};
				\addlegendentry{心形曲线}
			\end{axis}
		\end{tikzpicture}
		\caption{拟合样条效果}\label{ProblemE}
	\end{figure}
	

	
\end{document}